www.gusucode.com > 溷沌分析工具箱 - OpenTSTOOL1.11 > 混沌分析工具箱 - OpenTSTOOL1.11\mex-dev\Utils\mixembed.cpp
/***********************************************************************\
**
** Title: mixed state embedding
**
** File: mixembed.cpp
**
** Contents: main program
**
** Author: Kevin Bube - Drittes Physikal. Insitut Uni Goettingen
**
** Version 0.1: 4.10.2001
** 0.11: 5.10.2001 introduced type vector
** 0.3: 12.11.2001 added more embedding options
** 0.5: 24.5.2002 integrated into TSTOOL-CVS
**
** Bugs: Under-/Overflow checks are missing
**
** This is the implementation of a matlab mex-file. Its purpose is
** to embed time series into a mixed state vector.
** The invokation in matlab is mixembedalg (X, Y). For details type
** "help mixembed" at matlab prompt.
**
** Copyright 1997-2002 DPI Goettingen
** License http://www.physik3.gwdg.de/tstool/gpl.txt
**
\***********************************************************************/
#include "mex.h"
#include "matrix.h"
typedef struct{
double* data;
int rows, cols;
}matrix;
/*------------------------------------------------------------------------
*
* internal function 'embed'
*
* This function embeds time series into a mixed state vector.
*
* parameter: ts - mxn matrix of n time series of length m
* dims - 1xn matrix of embedding dimensions. One for each
* time series
* lags - 1xn matrix of time lags (in samples)
* shifts - 1xn-1 matrix which holds the amount of samples
* which the time series is shifted relatively to the
* first one
* startindex - pointer to an integer variable which holds the
* index of the time series value, which is the
* first value of the first embedding vector
* (this is used mainly for debugging purposes)
*
* return value: pointer to Matlab array of embedding vector
*
* preconditions: The dimensions of the time series must be the same. It
* Has to be checked if the values for the dimensions etc.
* are sane.
*
* side effects: -
*-----------------------------------------------------------------------*/
static mxArray *embed (matrix ts, matrix dims, matrix lags,
matrix shifts, int *startindex);
/*------------------------------------------------------------------------
*
* internal function 'max'
*
* This function returns the index of the maximum value of a vector. If
* there are more elements with the maximum value, the smallest index is
* returned.
*
* parameter: vector - pointer to double vector which should be searched
* length - integer value of the vector length
*
* return value: integer value of the index of the maximum element
*
* preconditions: -
*
* side effects: -
*-----------------------------------------------------------------------*/
static int max (const int *vector, int length);
/*------------------------------------------------------------------------
*
* internal function 'get_p_max_index'
*
* This function returns the first index of the first time series
* embedding vector.
*
* parameter: dims - 1xn dimension matrix
* lags - 1xn lags matrix
*
* return value: integer value of the maximum index
*
* preconditions: The arrays must be of the same lengths. The values of the
* dimensions must be nonegaive. The ones of the time lags
* must be positive.
*
* side effects: -
*-----------------------------------------------------------------------*/
static int get_p_max_index (matrix dims, matrix lags);
/**************************************************************************
*
* main function
*
*************************************************************************/
void mexFunction (int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
int i;
matrix ts, dims, lags, shifts;
int startindex;
int lags_need_free = false, shifts_need_free = false;
/* shifts needs init
*/
shifts.data = NULL;
shifts.rows = 0;
shifts.cols = 0;
/*-----------------------------------------------------------------------
* some simple sanity checks and data formating
*---------------------------------------------------------------------*/
switch (nrhs){
case 4:
if (!(mxIsDouble (prhs[3]) &&
!mxIsComplex (prhs[3]))){
mexErrMsgTxt ("corrupt input arg 4");
}
shifts.data = mxGetPr (prhs[3]);
shifts.rows = mxGetM (prhs[3]);
shifts.cols = mxGetN (prhs[3]);
/* fall through here */
case 3:
if (!(mxIsDouble (prhs[2]) &&
!mxIsComplex (prhs[2]))){
mexErrMsgTxt ("corrupt input arg 3");
}
lags.data = mxGetPr (prhs[2]);
lags.rows = mxGetM (prhs[2]);
lags.cols = mxGetN (prhs[2]);
/* fall through here */
case 2:
if (!(mxIsDouble (prhs[1]) &&
!mxIsComplex (prhs[1]) &&
mxIsDouble (prhs[0]) &&
!mxIsComplex (prhs[0]))){
mexErrMsgTxt ("corrupt input arg 1 or 2");
}
dims.data = mxGetPr (prhs[1]);
dims.rows = mxGetM (prhs[1]);
dims.cols = mxGetN (prhs[1]);
ts.data = mxGetPr (prhs[0]);
ts.rows = mxGetM (prhs[0]);
ts.cols = mxGetN (prhs[0]);
break;
default:
mexErrMsgTxt ("wrong number of input arguments");
}
/*
* set default values if required
*/
switch (nrhs){
case 4:
/* nothing to be done */
break;
case 3:
if (ts.cols > 1){
/* we don't need this stuff in case we have only one time series
*/
shifts.data = (double*) mxMalloc ((ts.cols-1) * sizeof (double));
if (shifts.data == NULL)
mexErrMsgTxt ("mexFunction(): error while allocating memory");
for (i = 0; i < ts.cols-1; i++)
shifts.data[i] = 0.0;
shifts.rows = 1;
shifts.cols = ts.cols-1;
shifts_need_free = true;
}
break;
case 2:
if (ts.cols > 1){
/* this would lead to an error if we had only one time series
*/
shifts.data = (double*) mxMalloc ((ts.cols-1) * sizeof (double));
if (shifts.data == NULL)
mexErrMsgTxt ("mexFunction(): error while allocating memory");
for (i = 0; i < ts.cols-1; i++)
shifts.data[i] = 0.0;
shifts.rows = 1;
shifts.cols = ts.cols-1;
shifts_need_free = true;
}
lags.data = (double*) mxMalloc (ts.cols * sizeof (double));
if (lags.data == NULL)
mexErrMsgTxt ("mexFunction(): error while allocating memory");
for (i = 0; i < ts.cols; i++)
lags.data[i] = 1.0;
lags.rows = 1;
lags.cols = ts.cols;
lags_need_free = true;
break;
default:
/* this should not happen */
mexErrMsgTxt ("internal error! please call the men in black");
}
/*
* check dimensions of the matrices we were given
*/
if (lags.rows != 1 ||
lags.cols != ts.cols ||
dims.rows != 1 ||
dims.cols != ts.cols){
mexErrMsgTxt ("wrong matrix dimensions");
}
if (ts.cols > 1){
/* skip this if we have only one time series
*/
if (shifts.rows != 1 ||
shifts.cols != ts.cols-1){
mexErrMsgTxt ("wrong matrix dimensions");
}
}
/*----------------------------------------------------------------------
* do the deed and format output
*--------------------------------------------------------------------*/
plhs[0] = embed (ts, dims, lags, shifts, &startindex);
plhs[1] = mxCreateDoubleMatrix (1, 1, mxREAL);
(mxGetPr (plhs[1]))[0] = (double) startindex;
/*
* clean up
*/
if (lags_need_free)
mxFree (lags.data);
if (shifts_need_free)
mxFree (shifts.data);
}
static mxArray *embed (matrix ts, matrix dims, matrix lags,
matrix shifts, int *startindex)
{
int i, j, k=0, l=0;
matrix p; /* embedding vector */
mxArray *mx_p;
int xindex;
/* get first index of time that should be predicted */
xindex = get_p_max_index (dims, lags);
/* determine size of Matlab array and create it */
/* We simply take the biggest shift when calculating p.rows.
This is not the optimal way as some time steps may be dropped
but in long time series this should not matter. */
p.rows = ts.rows - xindex - max ((int *)shifts.data, shifts.cols);
for (i = 0, p.cols = 0; i < ts.cols; i++)
p.cols += ((int) dims.data[i]);
mx_p = mxCreateDoubleMatrix (p.rows, p.cols, mxREAL);
if (mx_p == NULL)
mexErrMsgTxt ("embed(): error while allocating memory");
p.data = mxGetPr (mx_p);
/* finally it is time to get our hands dirty and fill the bucket */
for (i = 0; i < ts.cols; i++){
for (j = 0; j < dims.data[i]; j++){
/* there is no shift for the first time series */
if (i == 0){
for (k = xindex - (lags.data[i]*j) - 1;
k < ts.rows - lags.data[i]*j - 1;
k++){
p.data[l] = ts.data[k + i*ts.rows];
l++;
}
}
else{
for (k = xindex - (lags.data[i]*j) - 1 + shifts.data[i-1];
k < ts.rows - lags.data[i]*j - 1 + shifts.data[i-1];
k++){
p.data[l] = ts.data[k + i*ts.rows];
l++;
}
}
}
}
*startindex = xindex;
return mx_p;
}
static int get_p_max_index (matrix dims, matrix lags)
{
int *prod;
int i;
/* determine partial max. indizes */
prod = (int*) mxMalloc (dims.cols * sizeof (int));
if (prod == NULL)
mexErrMsgTxt ("get_p_max_index(): error while allocating memory");
for (i = 0; i < dims.cols; i++)
prod[i] = (dims.data[i] - 1) * lags.data[i] + 1;
i = max (prod, dims.cols);
i = prod[i];
mxFree (prod);
return i;
}
static int max (const int *vector, int length)
{
int index_max=0, i;
for (i = 0; i < length; i++){
if (vector[i] > vector[index_max])
index_max = i;
}
return index_max;
}